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1  Introduction 


The  objective  of  the  project  was  to  provide  means  for  establishing  reliable  quantitative  failure 
initiation  criteria  for  multiple  chip  modules  (MCMs),  electronic  packages,  laminated  composites 
and  adhesively  bonded  joints.  During  Phase  I,  an  easy-to-use,  reliable  and  robust  software  was 
developed,  with  a  graphic  user  interface,  based  on  the  innovative  methods  presented  in  [1],  [2] 
and  [3],  for  the  computation  of  generalized  flux/stress/thermal  intensity  factors  (GFIFs/GSIFs/ 
TSIFs)  and  the  strength  of  the  singularities  for  any  multi-material  interface  problem  involving  iso¬ 
tropic  or  anisotropic  materials,  subject  either  to  mechanical  or  thermal  loading,  in  a  two-dimen¬ 
sional  setting.  The  existing  software  product  Stress  Check  provided  the  framework  for  this 
development.  Stress  Check  is  based  on  the  p-  and  hp-version  of  the  finite  element  method,  capable 
of  a-posteriori  error  estimation  in  terms  of  the  data  of  interest. 

The  specific  accomplishments  during  the  Phase  I  project  are  summarized  below: 

•  Incorporation  of  the  modified  Steklov  formulation  presented  in  [1],  for  computing  the  strength 
of  the  singularities  (and  the  associated  eigenfunctions)  for  the  heat-transfer  and  elasticity  prob¬ 
lems,  into  the  software  product  Stress  Check. 

•  Numerical  solution  of  several  representative  test  cases  and  comparison  with  known  exact  solu¬ 
tions  to  demonstrate  the  robustness  and  accuracy  of  the  computation  of  eigenpairs. 

•  Implementation  of  the  algorithm  for  extracting  GFIFs  and  GSIFs.  This  algorithm,  based  on  the 
complementary  energy  principle  in  conjunction  with  the  p-  and  hp-versions  of  the  finite  ele¬ 
ment  method,  is  outlined  in  [2]. 
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•  Verification  of  the  accuracy  and  robustness  of  the  implemented  algorithm  by  computing  GFIFs 
and  GSIFs  for  problems  where  the  exact  solution  is  known. 

•  Implementation  of  an  algorithm  for  computing  the  thermal  stress  intensity  factors  (TSIFs), 
based  on  the  strategy  presented  in  [3]. 

•  Numerical  verification  of  some  thermo-elastic  problems  for  which  reference  solutions  are 
available  in  the  literature. 

•  Development  of  a  detailed  plan  for  the  implementation  of  this  technology  in  three  dimensions. 
The  three-  dimensional  setting  involves  vertex  singularities,  edge  singularities  and  vertex-edge 
singularities. 

The  successful  completion  of  these  activities  made  it  possible  to  address  the  following  impor¬ 
tant  questions: 

•  Which  are  the  characterizing  parameters  of  a  thermal,  elastic  and  thermo-elastic  solution  near  a 
singular  point  associated  with  multi-material  anisotropic  interface  problem  for  any  type  of  sin¬ 
gularity?  Consequently,  what  should  be  changed  in  the  design  so  as  to  minimize  the  likelihood 
of  failure  initiation?  Given  a  set  of  alternatives,  which  combination  of  materials  is  optimal,  and 
which  is  the  best  geometric  configuration? 

•  When  failure  is  observed  in  a  device,  how  should  it  be  modified  so  as  to  reduce  the  likelihood 
of  a  future  failure? 

•  How  accurate  are  the  parameters  that  influence  failure  initiation  obtained  by  the  numerical 
algorithm,  i.e.,  how  reliable  are  the  results  obtained  by  numerical  simulation?  How  does  the 
temperature  field  affect  the  solution  in  the  vicinity  of  the  singular  points? 

•  How  should  the  analyst  interpret  the  new  GFIFs/GSIFs  and  the  “eigenpairs”  when  correlating 
with  experimental  observations? 

•  How  to  extend  the  current  two-dimensional  capability  to  three  dimensions  where  the  full 
potential  of  the  technology  will  be  realized? 

A  detailed  description  of  each  activity  and  the  corresponding  formulations  are  given  in  the  next 

section. 

2  Accomplishments 


This  section  describes  the  activities  performed  during  the  Phase  I  project,  in  which  an  easy-to- 
use,  reliable  and  robust  software  was  developed  for  the  computation  of  the  generalized  flux/stress 
intensity  factors  and  the  strength  of  singularities  for  multi-material  interface  problems  subjected 
to  thermo-mechanical  loads  in  two-dimensions. 

The  solutions  of  linear  elastostatic  and  steady- state  heat  transfer  problems  in  the  vicinity  of 
crack  tips  were  an  intensive  subject  of  research  during  the  last  30  years.  Although  an  exact  solu- 
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tion  can  be  obtained  for  cracks  in  bodies  of  simple  geometries,  for  most  cases  involving  complex 
geometries,  anisotropic  materials,  and  cracks  at  bi-material  interfaces  only  a  numerical  approxi¬ 
mation  can  be  obtained.  Some  typical  singular  points  in  an  electronic  device,  for  example,  where 
failure  initiation  commonly  occurs,  are  illustrated  in  Figure  1 . 


FIGURE  1 .  Typical  sites  of  failure  initiation  in  an  electronic  device. 


The  solution  in  the  vicinity  of  singular  points  are  of  considerable  engineering  interest  (espe¬ 
cially  for  general  domains  containing  multi-material  interfaces,  and  anisotropic  materials) 
because  it  is  directly  or  indirectly  related  to  failure  initiation  in  composite  materials  and  electronic 
devices.  In  the  neighborhood  of  singular  points  the  exact  solution  of  two-dimensional  elastostatic 
problems  can  be  expanded  in  the  form: 


oo  OC. 

{uex)  =  A.  r  '{<!>. (6)}  (1) 

i  =  1 

where  {u^)  is  the  displacement  vector  in  the  x  and  y  directions,  r  and  0  are  polar  coordinates  cen¬ 
tered  on  the  singular  point;  ttj  are  called  eigenvalues  and  (t);(0)  are  called  eigenfunctions.  These 
eigenpairs  (ttj,  (]),•)  depend  on  the  material  properties,  the  geometry  and  the  boundary  conditions 
((|)j(0)  are  smooth  vector  functions).  The  A,-  are  coefficients  which  depend  on  the  loading.  Because 
of  their  close  analogy  to  stress  intensity  factors  in  linear  elastic  fracture  mechanics.  A,-  are  called 
generalized  stress  intensity  factors  (GSIFs).  In  the  case  of  linear  steady-state  heat  transfer  prob¬ 
lems,  the  solution  in  the  neighborhood  of  singular  points  is  analogous  to  Eq.  (1),  the  differences 
are  that  the  equation  is  in  a  scalar  form  and  the  coefficients  are  called  generalized  flux  intensity 
factors  ( GFIFs). 
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For  general  singular  points  the  exact  solution  {mjpj}  is  generally  not  know  explicitly,  i.e.,  nei¬ 
ther  the  exact  eigenpairs  nor  the  exact  GFEFs/GSIFs  are  known,  therefore  a  numerical  approxima¬ 
tion  must  be  found. 

The  stresses  in  the  same  neighborhood  can  be  computed  from  the  displacements  given  by  Eq. 
(1)  and  the  material  properties  as: 


oo  a.- 1 

{T)  =  £  A,.r  {^,-(6))  PI 

i  =  1 

where  V)/,(0)  depend  on  the  eigenfunctions  in  Eq.  (1)  and  the  material  coefficients.  It  is  clear  from 
Eq.  (2)  that  when  a,-  <  1,  the  stresses  become  singular  for  r=0. 

The  key  to  successful  failure  analysis  in  the  presence  of  singular  points  is  to  compute  reliably 
both  the  eigenpairs  and  the  GSIFs.  The  eigenvalues  characterize  the  strength  of  the  singularity,  the 
eigenfunctions  characterize  the  straining  modes  and  their  amplitudes  (the  GSIFs/GFIFs)  quantify 
the  amount  of  energy  residing  in  particular  straining  modes. 

A  general  method  for  computing  the  solution  in  the  vicinity  of  any  singular  point  has  been 
implemented  which  first  determines  the  eigenpairs,  followed  by  the  computation  of  the  GSIFs/ 
GFIFs  for  two-dimensional  problems.  During  the  Phase  I  project,  the  following  specific  objec¬ 
tives  were  achieved: 

1.  Incorporation  of  the  modified  Steklov  formulation  presented  in  [1],  for  computing  the 
strength  of  the  singularities  (and  the  associated  eigenfunctions)  for  the  heat-transfer  and 
elasticity  problems  in  two-dimensions,  into  the  software  product  Stress  Check. 

2.  Implementation  of  the  algorithm  for  extracting  the  GFIFs  and  GSIFs.  This  algorithm,  based 
on  the  Complementary  Energy  Principle  in  conjunction  with  the  p-  and  hp-  versions  of  the 
finite  element  method,  is  outlined  in  [2]. 

3.  Several  representative  test  cases  were  solved  to  demonstrate  the  robustness  and  accuracy  of 
the  computation  of  eigenpairs,  GFIFs  and  GSIFs  for  problems  where  the  exact  solution  was 
available. 

4.  Implementation  of  an  algorithm  for  computing  the  thermal  stress  intensity  factors  (TSIFs), 
based  on  [3]  and  performance  of  a  thermo-elastic  analysis  of  a  typical  problem. 

2.1  Computation  of  Eigenpairs 

The  implementation  of  the  modified  Steklov  method  into  Stress  Check  allows  for  the  compu¬ 
tation  of  the  eigenvalues,  and  the  corresponding  eigenvectors,  for  singularities  in  two-dimensional 
elastostatic  and  heat- transfer  problems.  Once  the  data  are  entered  into  the  program,  the  area 
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around  the  selected  singular  point  (but  excluding  the  singular  point)  is  internally  divided  into 
finite  elements  through  a  meshing  process  which  does  not  require  user  intervention.  This  ‘internal 
mesh’  is  arranged  in  a  circular  ring  around  the  singular  point  selected  by  the  user  in  such  a  way 
that  the  element  boundaries  coincide  with  the  material  interfaces.  The  number  of  elements  of  the 
internal  mesh  is  controlled  by  the  number  of  material  interfaces  around  the  singular  point  and  by 
the  solid  angle  of  the  singularity  as  shown  in  Figure  2.  The  largest  solid  angle  for  a  single  element 


FIGURE  2.  Typical  ‘internal  mesh’  around  singular  point  in  2D. 


is  limited  to  120°.  For  each  element  of  the  ‘internal  mesh’  we  compute  the  corresponding  stiffness 
and  mass  matrices.  Once  the  elemental  matrices  are  assembled  and  the  static  condensation  is  per¬ 
formed,  the  following  eigenvalue  problem  is  obtained: 

=  a[M]{««}  (3) 

where  [K^J  is  the  condensed  stiffness  matrix;  [M]  is  the  mass  matrix,  and  {uj^j  is  the  vector  of 
coefficients  corresponding  to  the  degrees  of  freedom  associated  with  the  circular  boundaries  of 
the  ‘internal  mesh’.  The  solution  of  the  eigenproblem  given  by  the  above  equation  yields  approxi¬ 
mation  for  the  eigenvalues  ai  and  the  corresponding  eigenvectors.  The  system  of  equations  is 
solved  for  increasing  polynomial  order  to  get  a  converging  sequence  of  solutions  (eigenpairs) 
using  routines  from  the  LAPACK  library. 

The  steps  and  fundamentals  for  obtaining  the  system  described  by  Eq.  (3)  are  described  in  the 
following. 
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Consider  a  domain  shown  in  Figure  3,  in  which  r,  0  are  the  coordinates  of  a  cylindrical 


FIGURE  3.  Solution  domain  and  notation  for  the  modified  Steklov  formulation. 


coordinate  system  located  in  the  singular  point.  By  formulating  the  weak  form  over  the  sin¬ 
gular  point  is  excluded  from  the  domain  of  interest  such  that  the  accuracy  of  the  finite  element 
solution  does  not  deteriorate  in  its  vicinity.  On  the  boundaries  and  r2  consider  either  traction- 
free  or  zero  displacement  boundary  conditions: 


-> 

T  =  0 


u  =  0 


on  i  =  1,  2. 


In  and  Uy  may  be  represented  as  follows: 


Using  Eq.  (5),  on  r3: 


l-M  =  ia/R)tc  (6) 

dr 

and  a  similar  condition  on  r4. 

Multiplying  the  equilibrium  equation  by  v  =  x  integrating 

using  Green’s  theorem,  and  following  the  steps  presented  in  [1],  the  modified  Steklov  weak  form 
is  obtained: 


Seek 


ae  C, 


O^ue  h\q*^)  X 
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B(m,v)  -  {Nji{u,v)  +  N^.{u,v))  =  a{Mif{u,v)  +  M^.iu,v)), 


V  ve 


(7) 


where  C  is  the  complex  plane  and 


B(m,v)  =  Jj([D]vH£]([Z)]M)^?Q, 

A/^(m,v)  -  |[Me][£]mM]r  =  /;  dQ, 
e 

N/uh  =  J[M(2][£]([£>^®^]«)]r=« 

0 

and  [D],  [Q]  and  [/?]  are  given  as  follows: 


1 

o 

(-sme)A 

0 

[D]  = 

dy 

0 

A  A 

_ 

(-sine)A^ 

r 

1 

COS0 

0 

COS0 

0  sin0 

[R]  = 

0 

sin0 

_  0 

sin0  COS0 

sin0 

COS0 

(8) 

(9) 

(10) 


(11) 


(12) 


and  [E]  is  the  material  matrix. 

Remark  1.  The  domain  does  not  include  singular  edge,  hence  no  special  refinement  of  the 
finite  element  mesh  is  required. 

Remark  2.  The  formulation  of  the  weak  form  was  not  based  on  the  assumption  that  the  material  is 
isotropic,  and  in  fact  can  be  applied  to  multi-material  anisotropic  interface. 

The  domain  is  divided  into  finite  elements  through  a  meshing  process,  as  described 
before.  The  polynomial  basis  and  trial  functions,  {<I>j),  are  defined  on  a  standard  element  in  the 
r\  space  such  that  -1<  ^  <1,  -1  <  ri  <1.  The  entries  of  the  unconstrained  stiffness  matrix  corre¬ 
sponding  to  b(u,v)  are  given  by  (see  Ref.  [1]): 

K,j  =  J  J  ([£>]{0,.})’'[£][D]{fl>.}rf£2  (13) 
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For  simplicity  traction-free  boundary  conditions  are  assumed  on  Fj  and  r2.  Considering  first 
iV(M,v)  =  N^iuh  +  N^.iu,v) ,  the  entries  of  the  matrix  [Nj^]  corresponding  to  the  bilinear  form  Nj^  are 
computed  using  Gauss  quadrature: 


5  3 

s  =  1  l,k  =  i 

where  and  are  the  weights  and  abscissas  of  the  Gauss  quadrature  points,  respectively,  and 
Pii  and  BP^^are  matrices  given  explicitly  in  [1]. 

The  entries  of  [M/j]  are: 


(0-0  3 

(My)«  =  X  PniWkPl^jiis)  (15) 

J  =  1  /,  i  =  1 

Expressions  similar  to  Eq.  (14)  and  Eq.  (15)  exist  for  the  matrices  and 

Denoting  the  set  of  all  coefficients  by  and  the  set  of  coefficients  associated  with  r3  and 

r4  by  {m^},  the  following  eigenproblem  is  obtained: 


{[K]  -  [yv„]  -  [)V^.]){«,.,}  =  a{[M„]  +  [M^.]){u^}  =  a[M]{u^}  (16) 

The  vector  which  represents  the  total  number  of  nodal  values  in  can  be  divided  into  two 
vectors  such  that  one  contains  the  coefficients  {w^},  the  other  contains  the  remaining  coefficients: 
{Ufot}^=  {  By  eliminating  {«(•„},  the  reduced  eigenproblem  is  obtained: 

=  a[M]{u^}  (17) 

Solution  of  the  eigenproblem  given  by  Eq.  (17)  yields  approximations  for  eigenpairs  with  high 
accuracy,  efficiency  and  robustness. 

A  new  module  inside  Stress  Check  allows  users  to  compute  the  eigenvalues  and  associated 
eigenfunctions  in  a  very  convenient  and  easy  to  use  way.  The  steps  necessary  to  compute  the 
eigenpairs  can  be  summarized  as  follows: 

•  First,  the  model  problem  is  loaded  into  Stress  Check  from  an  existing  neutral  file  or  it  is  created 
inside  the  program  using  the  existing  pre-processing  tools.  Stress  Check  is  a  p-version  finite 
element  analysis  program  developed  by  ESRD,  Inc.  for  the  solution  of  elastostatic  and  heat- 
transfer  problems.  It  has  a  Motif-based  graphic  user  interface  that  allows  for  the  creation  of  all 
the  geometric  features  of  the  problem  under  consideration.  After  the  geometry  was  created,  the 
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model  can  then  be  meshed,  the  material  properties  specified  and  the  boundary  conditions 
imposed.  Figure  4a  shows  a  typical  two-dimensional  model  problem  displayed  in  the  main 
window  of  Stress  Check. 


Display:  O  Report  if  Graph  O  File 


(^tion: 

None 

J-Integral 

Stress  Intensity  Factors 
Generalized  SIF 


If  Radius: 


#  of  terms : 


(a)  Main  Window  (b)  Fracture  Mechanics  diaiog  box 

FiGURE  4.  Stress  Check  Interface. 

♦  Second,  the  Fracture  Mechanics  dialog  box  is  loaded  by  a  menu  selection  from  the  main  win¬ 
dow  bar.  As  shown  in  Figure  4b,  the  fracture  mechanics  options  for  two-dimensional  elasticity 
are  the  J-integral,  the  Stress  Intensity  Factors  and  the  Generalized  Stress  Intensity  Factor.  The 
first  two  options  are  applicable  only  for  cracks  in  homogeneous  materials,  and  are  standard  fea¬ 
tures  of  Stress  Check.  The  Generalized  SIF  option  has  been  incorporated  as  part  of  this  Phase  I 
project 

•  Finally,  to  compute  the  eigenvalues  and  eigenvectors,  the  user  enters  the  number  of  eigenpairs 
to  be  computed  in  the  “#  of  terms”  region  of  the  dialog  box  and  then  selects  the  singular  point 
by  pointing  to  it  with  the  cursor  and  clicking  the  left  button  of  the  mouse.  The  results  are  dis¬ 
played  in  tabular  and/or  graphical  format  depending  on  the  display  selection  in  the  dialog  box. 

The  implementation  of  the  modified  Steklov  method  in  Stress  Check  was  tested  by  solving  a 
set  of  representative  benchmark  problems  for  which  exact  solutions  are  available  in  the  literature. 
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Two  groups  of  problems  were  investigated:  Steady  state  linear  heat-transfer  problems  and  linear 
elastostatic  2D-problems  as  described  in  the  next  sections. 

2.1 .1  Heat  transfer  problems 

Steady  state  linear  heat-transfer  problems  (also  called  scalar  problems)  in  the  neighborhood  of 
singular  points  are  considered  in  this  Section.  Here  u  denotes  the  temperature  field  in  a  domain. 
The  governing  differential  equation  is: 


a 


2  2  2 
d  u  ^  d  u  d  u 


=  0 


(18) 


where  aij  are  the  coefficients  of  heat  conduction  in  each  subdomain,  with  ciij=-ciji  and  the  satis¬ 
fying  the  elliptic  restriction:  aiin22  ’  <^12^  >  0  in  each  subdomain.  In  the  case  of  multimaterial 
interfaces,  it  is  assumed  that  the  materials  are  perfectly  bonded. 

Scalar  problem  1:  Isotropic  clamped-free  crack.  Circular  domain  of  unit  radius  with  a  crack 
along  the  positive  x-axis.  One  face  of  the  crack  (Tj)  has  zero  temperature  boundary  condition  and 
the  other  face  (r2)  is  flux  free.  The  outside  boundary  (Tr)  has  an  imposed  flux  (Figure  5) 


Boundary  Conditions: 

M  =  0  on  Tj ;  =  0  on  T; 

duldx  =  y  on  Tr 


FIGURES.  Scalar  problem  1.  Notation. 


The  exact  solution  for  this  problem  is  given  by  (ref.  [4]): 

«(r,e)  =  -1.35812/'"' sin(e/4)  +  0.970087/^‘'sin((3e)/4)  +  ...  (19) 

The  exact  values  of  the  first  two  eigenvalues  for  this  problem  are:  aj^=l/4  and  a2=3/4. 
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Scalar  problem  2:  Anisotropic  reentrant  comer.  Heat  transfer  problem  in  an  anisotropic  material 
governed  by  Eq.  (18)  with  aii=4,  ai2=0,  022=^-  The  boundary  consists  of  a  90°  reentrant  comer 
generated  by  two  flux-free  edges  T  j  and  r2,  which  meet  at  the  origin  of  the  coordinate  system  and 
u(0,0)=0  (Figure  6).  The  circular  boundary  of  the  domain  (Fi^)  is  loaded  by  flux  boundary  condi- 


Boundary  Conditions: 
u(0,0)  =  0 
dw/dO  =  Oon  Fj,  F2 


FIGURE  6.  Scalar  problem  2.  Notation. 


tion  which  corresponds  to  the  first  symmetric  eigenfunction  of  the  asymptotic  expansion  of  u(x,y) 
about  the  reentrant  corner  (Ref.  [9]): 

q^  =  (aiicos^0  +  a22cos^e)|^  +  ^sin2e(a22-«u)0|^]  (20) 

where  the  solution  M(r,0)  can  be  written  in  the  following  form: 

00  _9«  - 

u  =  2  (l  +  3sin  0)  cos  —  atan(2tan0)  ^2i) 

n  =  1 

The  GFBF  is  arbitrarily  selected  to  be  Ai=l,  while  the  others  are  A^=0,  /=2,  3, ...  The  exact 
values  of  the  first  two  eigenvalues  are:  a^2l?>  and  a2=4/3. 

Scalar  problem  3:  Internal  interface  with  two  materials.  Two  materials  perfectly  bonded  along  a 
common  edge  satisfying  the  following  equation: 

=  0  in  Q-  (22) 

with  the  following  flux  conditions  along  the  external  boundary  (Figure  7): 


du 

dr 


on  Ti=  dQi 


i  =  1,2 


(23) 
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The  material  coefficients  are:  pi  =  \Q,  P2=l\  the  eigenvalues  are:  a j^=0. 731691779 
a2=l .268308221.  and: 

(0)  =  I  cos[(l-a)e]  +  cisin[(l-a)0]^ort  Fj 
^  [  CjCOs[(l -a)0]  +  C2C3sin[(l -a)0]  ^  on  r2 


1*2(0) 


cos[(l  +  o)0]-C3sin[(l  +  a)0] on  Fj 
CjCOs[(l  +  fl)0]-C2C3sin[(l  +  o)0]  on  F2 


(25) 


where  ci=6. 31818181818182,  C2=-2. 68181818181818,  03=0. 64757612580273,  and 
a=0.26830822130025.  The  exact  solution  for  this  problem  is  (ref.  [5]): 


Boundary  Conditions: 
duldt  =  /(r,6j  on  Tj,  r2 


FIGURE  7.  Scalar  problem  3.  Notation. 


M(r,0)  =  A,r“'lii(0)+A2r“^M0)  (26) 

where  =  A2  =  1. 

Results:  Figure  8  shows  the  finite  element  meshes  used  for  the  scalar  problems.  For  problem  1, 
two  layers  of  geometrically  graded  elements  (with  a  common  factor  of  0.15)  towards  the  singular 
point  were  used.  For  problem  2  no  geometrically  graded  meshes  were  used,  while  for  problem  3 
only  one  layer  of  elements  was  placed  around  the  singularity. 

The  results  of  the  computation  of  the  eigenvalues  are  summarized  in  Table  2.  For  each  problem 
the  exact  value  of  the  first  two  eigenvalues,  as  obtained  from  the  corresponding  references,  and 
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FIGURE  8.  Mesh  and  boundary  conditions  for  scalar  probiems. 


those  computed  numerically  with  the  modified  Steklov  method  implemented  in  Stress  Check,  are 
included.  As  the  results  indicate,  the  correlation  between  numerical  and  exact  values  is  excellent. 


TABLE  1.  First  and  second  eigenvalues  for  the  scalar  model  problems 


Scalar 

Problem 

First  Eigenvalue 

Second  Eigenvalue 

Exact 

Numerical 

Exact 

Numerical 

1 

1/4 

0.250000000 

3/4 

0.750000000 

2 

2/3 

0.666666676 

4/3 

1.333333308 

3 

0.731691779 

0.731691779 

1.268308221 

1.268308223 

3 
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Finally,  Figure  9  shows  the  first  eigenfunction  as  computed  in  Stress  Check  for  each  scalar  prob¬ 
lem. 
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FIGURE  9.  First  eigenfunction  for  the  scalar  problems. 

2.1 .2  Elastostatic  problems 

Linear  elastostatic  model  problems  in  the  neighborhood  of  singular  points  are  considered  in 
this  section.  Here  u  =  {u^  Uy}^  denotes  the  displacement  vector  in  x,  y-directions  and  (5-^,  Gy  and 
Xxy  are  the  stresses.  For  multi-material  interfaces  continuity  of  displacements  and  tractions  across 
boundary  interfaces  is  assumed. 

Elastostatic  problem  1:  Traction-free  Isotropic  L-shaped  domain.  L-Shaped  plane  elastic  body 
(Figure  10)  loaded  along  the  boundaries  by  the  Mode  1  and  Mode  2  stress  components  obtained 
from  the  asymptotic  expansion  of  the  displacement  field  about  the  vertex: 

=  Aittir”'  '/i(0) 'AI®) 

Oy  =  A^a^r^'  '^i(e)+A2a2r“-  '^2(6) 
t  =  A,air“'  \,(e) +A2a2r“'  ’/!2(0) 
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where  A 1  and  A2  are  constants  analogous  to  the  mode  1  and  mode  2  stress  intensity  factors  in  lin¬ 
ear  elastic  fracture  mechanics;  0C|=Q.5444837368  and  a2~Q  9085291 898  are  the  first  and  second 
eigenvalues;  and^-,  g,-,  hi,  i=l,  2  are  functions  of  6  given  in  ref.  [6]. 


FIGURE  10.  Elastostatic  problem  1.  Notation. 


Elastostatic  problem  2:  Traction-free  crack  at  an  isotropic  material  interface.  Bi-material  inter¬ 
face  composed  of  two  homogeneous  and  isotropic  materials  with  continuity  of  tractions  and  dis¬ 
placements  across  the  interface  (Figure  11).  The  body  is  in  a  state  of  plane  strain  and  it  is  loaded 


FIGURE  11.  Elastostatic  problem  2.  Notation. 
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by  the  stress  field  corresponding  to  the  exact  solution  of  the  asymptotic  expansion  about  the  sin¬ 
gular  point; 


Tni^nr: 


1.5 


1  91  3  3 

— ={A^r[cos(£lnr)a^^+  sin(elnr)0^^]  +  /!r;;[cos(elnr)a^^+  (-sin 


(elnr))afj}L 


=  L5 


^^’6)1=  1,5 


-^r{A:;[cos(elnr)afe+  sin(elnr)afe]  + /s^//[cos(elnr)afe  + (-sin(elnr))afe]}| 
V37t 


where  and  a^Q  are  given  in  ref.  [7],  Ki  and  Kjj  are  the  stress  intensity  factors  and  e  is  given  by; 


e  = 


(3-4Vi)G2  +  Gi 

(3  -4v2)Gj  +  Gj 


=  0.07581178 


where  v  is  the  Poisson’s  ratio  and  G  is  the  shear  modulus  of  the  material.  The  first  two  eigenvalues 
for  this  problem  are  complex;  aj^  =  0.5  -i-  i  e  and  =  0.5  -  i  £. 

Elastostatic  problem  3:  Inclusion  problem.  Composite  body  consisting  of  two  dissimilar  isotro¬ 
pic,  homogeneous  and  elastic  wedges,  perfectly  bonded  along  the  interfaces  (Figure  12).  The 
body  is  loaded  by  the  stress  field  corresponding  to  the  exact  solution  of  the  asymptotic  expansion 
about  the  singular  point  as  given  in  ref.  [8].  The  eigenvalues  characterizing  the  stress  singularity 


FIGURE  12.  Elastostatic  problem  3.  Notation. 


at  (0,0)  are;  «!  =  0.512472160  and  =  0.730975740. 
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Results:  Figure  13  shows  the  mesh  used  in  solving  the  elastostatic  model  problems.  For  problems 
1  and  2,  two  layers  of  geometrically  graded  elements  (with  a  common  factor  of  0.15)  towards  the 
singular  point  were  used.  For  problem  3  only  one  layer  of  elements  was  placed  around  the  singu¬ 
larity. 


Problem  1 


FIGURE  13.  Mesh  and  boundary  conditions  for  elastostatic  problems. 
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The  results  of  the  computation  of  the  eigenvalues  are  summarized  in  Table  2.  For  each  prob¬ 
lem,  the  exact  value  of  the  first  two  eigenvalues  as  obtained  from  the  corresponding  references, 
and  those  computed  numerically  with  the  algorithm  implemented  in  Stress  Check  are  included.  As 
the  results  indicate,  the  correlation  between  numerical  and  exact  values  is  excellent. 


TABLE  2.  First  and  second  eigenvalues  for  the  elastostatic  model  problems 


Elastostatic 

Problem 

First  Eigenvalue 

Second  Eigenvalue 

Exact 

Numerical 

Exact 

Numerical 

1 

0.5444837368 

0.5444837375 

0.9085291898 

0.9085291893 

2 

0.5+0.07581178/ 

0.5+0.07581178/ 

0.5-0.07581178/ 

0.5-0.07581178/ 

3 

0.512472160 

0.5124721606 

0.730975740 

0.7309757404 

Finally,  Figure  14  shows  the  first  set  of  eigenfunctions  corresponding  to  problems  1  and  3. 
They  include  the  eigenfunctions  associated  with  the  two  displacement  components  Uy,  and  the 
eigenstresses  <5y  x^. 

2.2  Computation  of  Generalized  Flux/Stress  Intensity  Factors 

Once  the  eigenpairs  have  been  computed,  they  are  used  for  extracting  the  GFIFs/GSIFs  from 
the  finite  element  solution.  The  procedure  is  described  for  the  case  of  elasticity.  The  case  of  heat 
transfer  is  analogous. 

After  solving  the  elastostatic  problem  over  the  entire  domain  Q.  by  means  of  the  finite  element 
method  based  on  the  displacement  formulation,  ufe  is  obtained.  A  small  sub-domain  around  the 
singular  point  is  considered  next.  Defining  as  the  set  of  interior  points  of  a  circle  of  radius  R, 
centered  on  the  point  P,  is  defined  by  n  5/j,  and  Fr  is  the  circular  part  of  its  boundary.  See 
Figure  15. 

The  complementary  variational  principle  over  Qr  can  be  stated  as: 

Seek  Cq  g  such  that 

5,((To,o,)  =  F,(o,)  V  ai  e  (27) 
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Problem  3 


FIGURE  14.  Eigenfunctions  for  elastostatic  problems. 

E^{Q.j^)  being  the  statically  admissible  space  (see  a  detailed  definition  in  ref.  [2]),  and  and 
are  given  by: 


5,(ao,o,)  =  IJaoier'aiJQ  (28) 

F,(ai)  =  I  V[Q]aids  (29) 

an',"’ 

where  [£]  is  the  material  matrix,  is  that  part  of  the  boundary  where  the  displacement  vector 

u  is  prescribed,  and  [Q]  is  given  in  Eq.  (12). 
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For  the  complementary  weak  form  the  trial  and  test  spaces  are  chosen  to  be  linear  combina¬ 
tions  of  the  eigenstresses,  which  are  computed  from  the  eigenpairs,  using  the  stress-strain  rela¬ 
tionship  and  Hooke’s  law.  The  ij-th  term  of  the  compliance  matrix  which  corresponds  to  the 
bihnear  form  in  Eq.  (27)  is  given  by: 


(a.^  +  a^-l) 


{Du)kE^i{Du)idrd% 


being  the  real  part  of  the  eigenpair,  and  (Duf  is  a  vector  corresponding  to  the  j-th  eigenfunc¬ 
tion  and  is  given  in  ref.  [2]. 

The  eigenstress  tensor,  having  been  derived  from  the  eigenpairs,  automatically  satisfies  the 
boundary  conditions  on  all  boundaries  except  Tr,  so  that  the  linear  form  in  Eq.  (29)  degenerates 
into  an  integral  over  the  circular  boundary  Tr  alone. 

Replacing  the  vector  u  in  Eq.  (29)  with  the  approximated  finite  element  solution  up^  on  r3 
(see  Figure  3),  the  j-th  term  of  the  load  vector  corresponding  to  the  linear  form  Eq.  (29)  becomes: 

(FJ  =t?|  uo  [Q][E\[Du^^]  dQ  (31) 

•^01  ^  ^  {r^R) 

1  FE  FE  FE  FE  ^ 

where  Uo  =  COS0  +  W^  sin0),(Mj^  sin0  +  w^  cos0))  . 

Solving  Eq.  (27),  an  approximation  for  the  coefficients  of  the  asymptotic  expansion  (the 
GSEFs)  is  obtained.  Numerical  tests  demonstrated  that  the  rate  of  convergence  of  the  GSIFs  is  as 
fast  as  the  convergence  of  the  strain  energy,  therefore  the  method  is  “superconvergent”. 
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The  computation  of  the  GFIFs/GSIFs  was  integrated  in  Stress  Check  within  the  same  module 
described  in  Section  2.1  for  the  computation  of  eigenvalues  and  associated  eigenfunctions  for 
two-dimensional  problems  in  heat  transfer  and  elasticity.  The  steps  to  compute  the  GFIFs/GSIFs 
are  as  follows: 

•  The  model  problem  is  loaded  into  Stress  Check  from  an  existing  neutral  file,  or  it  is  created 
inside  the  program  using  the  existing  pre-processing  tools. 

•  After  executing  the  corresponding  analysis  (elasticity  or  heat  transfer,  depending  on  the  prob¬ 
lem),  the  Fracture  Mechanics  dialog  box  is  loaded  by  menu  selection  from  the  main  window 
bar. 

•  The  GFIFs/GSIFs  is  computed  by  entering  the  radius  of  the  extraction  circle,  the  number  of 
terms  to  be  computed  in  the  “#  of  terms”  region  of  the  dialog  box,  and  then  selecting  the  singu¬ 
lar  point  by  pointing  to  it  with  the  mouse  cursor.  The  results  are  displayed  in  tabular  and/or 
graphical  format,  depending  on  the  display  selection  in  the  dialog  box. 

A  set  of  benchmark  problems  for  which  analytical  (exact)  solutions  are  known  were  investi¬ 
gated.  These  problems  have  been  selected  as  representative  of  the  types  of  singularities  present  in 
practical  engineering  situations  associated  with  linear  steady  state  heat  transfer  and  elastostatic 
models.  The  six  problems  described  in  Section  2.1.1  and  Section  2.1.2  are  considered  for  the  com¬ 
putation  of  GFIFs  (scalar  problems)  and  GSIFs  (elastostatic  problems). 

TABLE  3.  First  and  second  GFIFs  for  the  scalar  problems 


Scalar 

Problem 

Ai 

A2 

Exact 

Numerical 

Exact 

Numerical 

1 

1.358097 

1.328489 

0.970087 

0.970085 

2 

1.0 

0.998462 

0.0 

0.000028 

3 

1.0 

1.001243 

1.0 

1.001222 

Table  3  shows  the  first  and  second  generalized  flux  intensity  factors  (Aj,  A2)  for  the  heat 
transfer  problems,  as  computed  using  the  present  algorithm  and  the  corresponding  exact  values. 
The  convergence  characteristics  of  the  extraction  procedure  are  illustrated  in  Figure  16,  which 
shows  the  values  of  Aj  and  A2  as  a  function  of  the  number  of  degrees  of  freedom  (DOF)  for  scalar 
problem  1.  As  the  number  of  degrees  of  freedom  of  the  finite  element  solution  is  increased,  the 
GFIFs  reach  a  limit  value  which  is  practically  independent  of  the  discretization. 


21 


Accomplishments 


Scalar  isotropic  clainped-rree  crack 
Solution  =  SOLj  runs  1  to  8  (Up=16) 
Generalized  Flux  Intensity  Factors;  Int.  Radius  =  0.9 


FIGURE  16.  Convergence  of  A1  and  A2  for  scalar  problem  1. 


Table  4  shows  the  first  and  second  generalized  stress  intensity  factors  (A|,  A2)  for  the  elasto- 
static  problems.  Also  included  are  the  corresponding  exact  values.  Figure  17  shows  the  conver- 

TABLE  4.  First  and  second  GSIFs  for  the  elastostatic  model  problems 


Elasticity 

Problem 

Ai 

A2 

Exact 

Numerical 

Exact 

Numerical 

1 

1.0 

0.999699 

1.0 

0.999989 

2 

1.0 

1.000122 

1.0 

0.999632 

3 

2.506628 

2.508634 

i 

2.506628 

2.506899 

gence  of  the  GSIFs  for  elasticity  problem  1  as  obtained  from  Stress  Check  in  tabular  and  graphical 
forms. 
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FIGURE  17.  Convergence  of  A1  and  A2  for  elasticity  problem  1. 


2.3  Computation  of  Thermal  Stress  Intensity  Factors 

When  an  elastic  body  is  subjected  to  thermal  loading,  the  first  two  coefficients  of  the  asymp¬ 
totic  expansion  about  the  singular  points  in  Eq.  (1)  are  called  the  thermal  stress  intensity  factors 
(TSIFs).  The  computation  of  the  TSIFs  can  be  performed  by  using  the  modified  Steklov  method 
for  the  computation  of  the  eigenpairs;  the  minimum  complementary  energy  principle  to  obtain  the 
stress  intensity  factors;  and  Richardson’s  extrapolation  to  determine  the  TSIFs  as  the  limiting  pro¬ 
cess  when  ^  0  (Figure  3).  The  details  of  the  formulation  are  available  in  [3], 

To  compute  the  TSIFs,  a  numerical  algorithm  was  integrated  into  Stress  Check  within  the 
same  interface  as  the  one  used  for  the  computation  of  the  GSIFs.  The  implementation  involves  the 
following  steps: 

Step  1.  Solve  the  thermal  problem:  Find  the  temperature  distribution  over  the  entire  domain  given 
the  flux  and  temperature  boundary  conditions.  Compute  the  smallest  eigenvalue  (Pj)  asso¬ 
ciated  with  the  flux  singularity  by  the  modified  Steklov  method  as  described  in  Section 
2.1. 
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step  2.  Solve  the  thermoelastic  problem:  Using  the  temperature  distribution  obtained  from  the 
solution  of  the  thermal  problem  as  input,  obtain  the  corresponding  displacement  field  of 
the  elasticity  problem. 

Step  3.  Compute  eigenvalues  and  SlFs:  Extract  the  first  two  eigenvalues  (ccj,  Oq)  and  the  first  two 
stress  intensity  factors  (Kj,  Kjj)  from  the  thermoelastic  problem  for  a  given  radius  (R)  of 
the  integration  circle.  Stress  Check  will  automatically  compute  the  stress  intensity  factors 
for  three  values  of  the  integration  circle:  0.9/?,  R  and  1.1/?. 

Step  4.  Compute  TSIFs:  Using  Richardson’s  extrapolation,  project  the  three  values  of  each  stress 
intensity  factor  computed  in  Step  3  to  R=0.  The  error  in  the  computation  of  the  stress 
intensity  factors  is  of  the  order  f/?)^,  where  q=^j-aj+l  for  the  first  intensity  factor  (Kj) 
and  ^=P;-a2+l  the  second  intensity  factor  (Kjj). 

The  implementation  is  illustrated  with  the  solution  of  a  representative  thermo-elastic  problem 
in  two-dimensions. 

Central  crack  in  a  rectangular  plate:  Consider  a  rectangular  isotropic  plate  subjected  to  two 
different  thermal  loadings  for  which  numerical  results  are  available  in  the  literature.  The  plate  of 
width  2W,  length  2L  and  crack  of  length  2a-2.0  with  L/W=1.0  and  a/W-0.2  is  solved  for  two 
thermal  loadings  representing  modes  I  and  II  respectively  (Figure  18).  Since  the  results  are  inde- 


A 


B 


ik 

Mode  I  Loading 

AB,  BC,  CD,  DA:  T=100 
EF:  T=0 

2L 

V 

Mode  II  Loading 

AB  :  T=100 

CD:  T=-100 

BC,  DA,  EF:  q„=0 

FIGURE  18.  Plate  with  central  crack.  Notation. 


pendent  of  the  thermal  conductivity  for  isotropic  materials,  the  heat  conduction  coefficients  are 
taken  as:  aij=a22=^-0,  aj2=0.  The  mechanical  properties  of  the  material  are:  Modulus  of  Elastic¬ 
ity  £=1.0;  Poisson’s  ratio  v=0.3;  coefficient  of  thermal  expansion  a=0.01.  Plane-strain  condi¬ 
tions  are  assumed.  Because  of  symmetry,  only  one  half  of  the  plate  is  analyzed.  Figure  19  shows 
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the  finite  element  mesh  used  to  solve  the  problem  for  mode  I  and  mode  n  loadings.  Only  one  layer 
of  geometrically  graded  elements  towards  the  crack  tip  was  used  in  this  case. 


The  values  of  the  stress  intensity  factor  Kj  ^^-^computed  for  three  radii  of  the  integration  circle 
for  Mode  I  loading  are  shown  in  Table  5,  together  with  the  values  computed  using  Richardson’s 

TABLE  5.  TSIFs  for  Mode  I  loading.  Computed  and  extrapolated  values. 


R 

K/l) 

k/2) 

0.45 

1.283830 

0.809215 

0.50 

1.336565 

0.811445 

0.79915 

0.55 

1.389077 

extrapolation  and  K^^).  For  Mode  I  loading  Kji=0.  The  smallest  eigenvalue  associated  with 
the  flux  singularity  was  computed  to  be  |37=0.5  and  the  first  two  eigenvalues  associated  with  the 
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elasticity  singularity  are:  ai=a2=0.5.  Figure  20  shows  the  convergence  of  the  mode  I  TSIF  as 
obtained  from  Stress  Check  in  tabular  and  graphical  forms,  as  a  function  of  the  number  of  degrees 
of  freedom  (DOF). 


Hxtr^olated  TSIF 

DOF  TSIFl 

TSIF2 

53 

5.815359e-01 

1.137471e-14 

155 

7.314151e-01 

9.745229e-16 

273 

7.524530e-01 

hl.761370e-15 

439 

7.785932e-01 

7.229132e-15 

653 

7.962741e--01 

1.814041e-14 

915 

7.994875e-01 

2.1992690-14 

1225 

7.988002e-01 

5.5380850-15 

1583 

7.991494e-01 

5.1265670-15 

FIGURE  20.  Convergence  of  Mode  I  TSIF. 


The  values  of  the  stress  intensity  factor  Ku  ^^-^computed  for  3  radii  of  the  integration  circle  for 
Mode  II  loading  are  shown  in  Table  6,  together  with  the  values  computed  using  Richardson’s 


TABLE  e.TSIFs  for  Mode  II  loading.  Computed  and  extrapolated  values. 


R 

Ku® 

K„® 

0.45 

0.027258 

0.122953 

0.50 

0.016626 

0.12170 

0.123231 

0.55 

0.005965 

extrapolation  {Ku  and  For  Mode  II  loading  Ki=0. 
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Finally,  the  extrapolated  values  of  the  thermal  stress  intensity  factors  Kj  (for  Mode  I  loading) 
and  Kjj  (for  Mode  II  loading)  are  compared  with  the  values  published  in  the  literature  in  Table  7. 

table7.TSIFs  for  Mode  I  and  II  loading.  Comparison  with  references. 


Present 

Ref  [10] 

Ref  [11] 

Ref  [12] 

Method 

Ki 

0.8593 

0.7759 

0.7759 

0.7992 

Kii 

0.1317 

0.1207 

0.1185 

0.1217 

The  numerical  example  for  a  crack-tip  singularity  indicates  that  the  method  implemented  in 
Stress  Check  for  the  computation  of  the  thermal  stress  intensity  factors  is  accurate  and  efficient. 

3  Extension  to  Three-dimensions 


The  proof  of  concept  completed  for  two-dimensional  problems  can  be  extended  to  three- 
dimensional  problems  that  typically  occur  in  the  design  and  manufacture  of  electronic  compo¬ 
nents.  This  involves  detailed  development  and  implementation  of  algorithms  for  the  computation 
of  the  eigenpairs  that  characterize  the  temperature  and  displacement  fields  in  the  vicinity  of  singu¬ 
larities  caused  by  multi-material  interfaces  and  certain  topological  details  in  three  dimensions; 
extraction  of  the  generalized  flux  and  stress  intensity  factors  in  three  dimensions,  and  a  posteriori 
estimation  of  the  error  in  the  computed  flux  and  stress  intensity  factors. 

3.1  Objectives 

The  main  goal  of  Phase  I  was  to  develop  easy-to-use,  reliable  and  robust  software,  with  a 
graphical  user  interface,  based  on  the  innovative  methods  presented  in  Section  2,  for  the  computa¬ 
tion  of  GFIFs/GSIFs,  TSIFs  and  the  strength  of  singularities  for  any  multi-material  interface  prob¬ 
lem  involving  isotropic  or  anisotropic  materials,  subjected  either  to  mechanical  or  thermal 
loading,  in  a  two-dimensional  setting.  The  existing  software  product  Stress  Check  provided  the 
framework  for  this  development.  The  following  specific  objectives  are  considered  as  an  extension 
of  the  methods  implemented  during  the  Phase  I  project: 

(1)  Develop  and  implement  an  algorithm  for  the  computation  of  eigenpairs  that  characterize 
the  temperature  distribution  in  the  vicinity  of  singular  vertices  and  edges  in  multi-material 
interface  problems:  Investigate  the  extension  of  the  Steklov  method  developed  for  two- 
dimensional  problems  for  the  computation  of  eigenpairs  in  three-dimensions.  Incorporate 
into  Stress  Check  the  modified  Steklov  formulation  for  computing  the  strength  of  singu¬ 
larities  (and  the  associated  eigenfunctions)  for  the  heat-transfer  problem. 
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(2)  Develop  and  implement  an  algorithm  for  the  computation  of  the  generalized  flux  intensity 
factors  (GFIFs)  in  the  case  of  steady  state  heat  conduction  problems:  The  algorithm  is 
based  on  the  Complementary  Energy  Principle  in  conjunction  with  the  p-  and  hp-versions 
of  the  finite  element  method. 

(3)  Develop  and  implement  an  algorithm  for  the  computation  of  eigenpairs  that  characterize 
the  strain  distribution  in  the  vicinity  of  singular  vertices  and  edges  in  multi-material  inter¬ 
face  problems:  The  modified  Steklov  formulation  for  computing  the  strength  of  singulari¬ 
ties  (and  the  associated  eigenfunctions)  can  be  extended  for  elasto-static  problems  in 
three-dimensions 

(4)  Develop  and  implement  an  algorithm  for  the  computation  of  the  generalized  stress  inten¬ 
sity  factors  (GSIFs)  when  an  elastic  body  is  subjected  to  mechanical  loads:  The  algorithm 
is  based  on  the  Complementary  Energy  Principle  in  conjunction  with  the  p-  and  hp-ver- 
sions  of  the  finite  element  method. 

(5)  Develop,  implement,  test  and  document  an  algorithm  for  computing  the  thermo-elastic 
stress  intensity  factors  (TSIFs)  in  three-dimensions. 

(6)  Develop  an  algorithm  for  the  coupled  thermo-elastic  problem  accounting  for  dependence 
of  material  properties  on  temperature  in  two-dimensions.  This  will  involve: 

-  Computation  of  the  nonlinear  heat  transfer  solution  and  investigation  of  the  decompo¬ 

sition  of  the  solution  into  asymptotic  terms. 

-  Extraction  of  the  generalized  flux  intensity  factors  and  generalized  stress  intensity 
factors  for  non-constant  (temperature  dependent)  material  properties. 

-  Test  implementation  and  documentation  of  the  algorithms. 

Successful  completion  of  these  activities  will  address  the  following  important  questions  for 

the  most  general  (three-dimensional)  problem: 

•  What  is  the  thermal,  elastic  and  thermo-elastic  solution  near  a  3D  singular  point  associated 
with  vertex,  edge,  and  vertex-edge  multi-material  anisotropic  interface  problem?  Conse¬ 
quently,  what  should  be  changed  in  the  design  so  as  to  minimize  the  likelihood  of  failure  initia¬ 
tion?  Given  a  set  of  alternatives,  which  combination  of  materials  is  optimal,  and  which  is  the 
best  geometric  configuration?  Which  singularity  is  worse  from  th  point  of  view  of  durability: 
The  vertex,  the  edge,  or  the  vertex-edge  singularity? 

•  When  failure  is  observed  in  a  device,  how  should  it  be  modified  so  as  to  reduce  the  likelihood 
of  a  future  failure?  Which  mode  of  failure  is  to  be  expected  next? 

•  How  accurate  are  the  parameters  that  influence  failure  initiation  obtained  by  the  numerical 
algorithm,  i.e.,  how  reliable  are  the  results  obtained  by  numerical  simulation? 

•  How  does  the  temperature  field  affect  the  solution  in  the  vicinity  of  singular  points  in  three- 
dimensions? 
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•  How  should  the  analyst  interpret  the  new  GFIFs/GSIFs  and  the  ‘  ‘eigenpairs”  when  correlating 
with  experimental  observations? 

•  How  do  the  temperature  dependent  materials  influence  the  singular  solution?  Consequently, 
what  should  be  done  to  prevent  or  reduce  the  likelihood  of  failure  initiation  and  propagation? 

Following  is  a  detailed  description  of  the  plan  for  the  development  and  implementation  of  the 
activities  described  above. 

3.2  Implementation  Plan 

Three-dimensional  singularities  are  considerably  more  difficult  to  analyze  than  two-dimen¬ 
sional  ones,  where  only  one  type  of  singularity  exists.  In  3D  in  a  neighborhood  of  the  edges  and 
the  vertices  the  singular  behavior  is  different. 

Edge  Singularities:  If  a  coordinate  system  (x,y,z)  is  located  at  an  edge,  with  the  y-axis  along 
the  edge,  then  there  are  three  edge  GSIFs  which  are  y-dependent:  A/yj,  Aji(y)  andAjj^y).  These 
edge  GSIFs  are  analytic  along  the  edge,  however  they  become  singular  when  this  edge  intersects 
with  a  free  plane,  at  a  vertex.  In  the  neighborhood  of  an  edge-vertex  type  geometry  the  GSIFs  can 
be  represented  once  again  by  vertex  and  vertex-edge  stress  intensity  factors.  For  example, 
A,(j)  =  +  smoother  terms. 

/ 

Vertex  Singularities:  In  the  neighborhood  of  a  vertex,  and  away  from  edge-vertex  geometry, 
the  displacement  field  can  be  represented  by  only  one  vertex  intensity  factor  and  the  correspond¬ 
ing  eigenpairs.  Investigating  the  mathematical  behavior  of  the  singularities  in  3-D  is  an  active 
field  of  research  in  the  mathematical  community,  and  the  decomposition  of  the  displacement  field 
into  singular  and  regular  terms  is  documented  in  some  recent  papers.  The  application  of  the  modi¬ 
fied  Steklov  method  and  the  complementary  energy  principle  for  extracting  GSIFs  in  3D  will  be 
addressed. 

3.2.1  Introduction  and  Notation 

The  solutions  of  linear  heat-transfer  problems  in  three-dimensions,  for  example  in  the  vicinity 
of  any  singular  point,  can  be  decomposed  into  three  different  forms,  depending  whether  the  singu¬ 
lar  point  is  in  the  neighborhood  of  an  edge,  a  vertex  or  an  intersection  of  the  edge  and  the  vertex. 
Mathematical  details  of  the  decomposition  can  be  found  e.g.  in  [13]-[16]  and  the  references 
therein.  A  representative  three-dimensional  domain  denoted  by  ^2,  which  contains  typical  3D  sin¬ 
gular  points  is  shown  in  Figure  21.  Vertex  singularities  arise  in  the  neighborhood  of  the  vertices 
A,-,  and  the  edge  singularities  arise  in  the  neighborhood  of  the  edge  singularities  Ay.  Close  to  the 
vertex/edge  intersection,  vertex-edge  singularities  arise.  Of  interest  is  the  solution  of  the  mathe¬ 
matical  problem: 


=  0  inQ 


(32) 
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Aj-  Vertex  / 

Ar  Edge  between  >4/ and  Aj 


FIGURE  21.  Typical  3D  singularities. 


M  =  0 

on  Fq  c  do. 

(33) 

0 

dn 

on 

(34) 

where  u(xj,  X2,  Xj)  denotes  the  temperature  field  (in  the  following  xj,  X2  and  xj  will  be  either  Car¬ 
tesian,  cylindrical  or  spherical  coordinates),  and  r^,  u  =  3Q .  It  shall  be  assumed  the  edges  that 
intersect  at  vertices  are  not  curved,  and  that  crack  faces,  if  any,  lie  in  a  plane. 

Edge  Singularities:  The  edges  denoted  by  Ay,  which  connect  the  vertices  A,-  and  Aj,  are  consid¬ 
ered  first.  Moving  away  from  the  vertex  a  distance  6/2,  and  creating  a  cylindrical  subdomain  of 
radius  r  =  e  with  the  edge  Ay  as  its  axis,  a  subdomain  is  defined  in  the  vicinity  of  the  edge  denoted 
by  Eg  g(  Ay) .  Figure  22  shows  the  edge  singularity  subdomain  Eg  £(Ai2). 

The  solution  in  85  -  can  be  decomposed  as  follows: 

^  ^  »  (35) 

'■"‘(In'')  iA6)  +  v(r,  e,  z) 

k=ls  =  0 

where  5  >  0  is  an  integer  which  is  zero  unless  a/^  is  an  integer,  are  called  edge  eigenval¬ 

ues;  ai^(z)  are  called  the  edge  flux  intensity  functions  (EFIFs),  are  analytic  in  z  but  can  become 
very  large  as  they  approach  one  of  the  vertices;  and  fks(Q),  called  eigenfunctions,  are  analytic  in  0. 
The  function  v(r,  e,z)  belongs  to  H^(e).  It  is  assumed  that  for  k  <  is  not  an  integer,  therefore 
Eq.  (35)  becomes; 
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A2 


FIGURE  22.  The  edge  neighborhood  eg,e(A,2). 


K 

u{r,Q,z)  =  +  v(r,  0,  z).  (36) 

it=  1 

Vertex  Singularities:  A  ball  of  radius  p  =  5,  centered  in  the  vertex  Aj  for  example,  is  con¬ 
structed  and  intersected  by  the  domain  Q.  Then,  a  cone  having  an  opening  angle  (|)  =  a  is  con¬ 
structed  along  every  edge  intersecting  at  Aj,  and  removed  from  the  previously  constructed 
subdomain,  as  shown  in  Figure  23.  The  resulting  vertex  subdomain  is  denoted  by  VgfAjj,  and  the 
solution  u  can  be  decomposed  in  FgfAjj  using  a  spherical  coordinate  system  by: 

•>  ^  ®  ^  > 

M(p,(|),  0)  =  IS  b,q  p^'(lnp)%^((|),  0)  +  v(p,  (]),  0)  (37) 

l=lq=0 

where  g  ^  0  is  an  integer  which  is  zero  unless  Ji  is  an  integer,  >  Ji  are  called  vertex  eigenval¬ 
ues,  and  hiq(^,Q) ,  called  the  eigenfunctions,  are  analytic  in  (})  and  6.  The  big  are  called  vertex  flux 
intensity  factors  (VFIFs).  The  function  v(r,  0,z)  belongs  to  H^(V).  It  is  assumed  that  Ji  for  /  <  L  is 
not  an  integer,  therefore,  Eq.  (37)  becomes: 

«(p,^,0)  =  p^'  m,e)  +  v(p,^,0)  (38) 

/=  1 
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FIGURE  23.  The  vertex  neighborhood  VgCAi). 


Vertex-Edge  Singularities:  The  most  complicated  decomposition  of  the  solution  arises  in  case 
of  vertex-edge  intersections.  For  example,  consider  the  neighborhood  where  the  edge  A.\2 
approaches  the  vertex  Aj.  A  spherical  coordinate  system  is  located  in  the  vertex  A^,  and  a  cone 
having  an  opening  angle  (j)  =  a  with  its  vertex  coinciding  with  A  j  is  constructed  with  A12  being  its 
center  axis.  This  cone  is  terminated  by  a  ball-shaped  basis  having  a  radius  p  =  5,  as  shown  in 
Figure  24.  The  resulting  vertex-edge  subdomain  is  denoted  by  Ve5  £(Aj,  A12),  and  the  solution  u 
can  be  decomposed  in  F85  gCAj,  A12): 

K  S  f  L  ^  L  Q  ^  ^ 

M(p,(t),e)  =  P^'  +  f^ksiP)  (sin(t))“‘[ln(sin<D)]VAe)+^]^C/qP^'(lnp)’/j/9(<l),e)  +  v(p,<l),e) 


k=ls  =  0\J  =  \ 


/=  1?  =  0 


(39) 


where  mjf^(p)  is  analytic  in  p;  gks{Q)  is  analytic  in  0,  and  higi(^,B)  is  analytic  in  (j)  and  9.  The  func¬ 
tion  v(r,  e,z)  belongs  to  H^(Ve).  Again  it  is  assumed  that  Ji  for  /  <  L  is  not  an  integer,  and  for  k 
<K  is  not  an  integer,  therefore,  Eq.  (39)  becomes: 


m(p,(|),  0) 


k=  1 


/=  1 


^  L 

{sin^f'gkW  +  '^Ci  p 
/=  1 


hi(<^,Q)  +  v(p,(j),  9) 


(40) 


The  eigenvalues  and  the  eigenfunctions  are  associated  pairs  (eigenpairs)  which  depend  on  the 
material  properties,  the  geometry,  and  the  boundary  conditions  in  the  vicinity  of  the  singular  point/ 
edge  only.  Similarly,  the  solution  for  problems  in  linear  elasticity,  in  the  neighborhood  of  singular 
points/edges  is  analogous  to  Eq.  (35)-Eq.  (40),  the  differences  are  that  the  equations  are  in  vector 
form  and  the  eigenpairs  may  be  complex.  For  general  singular  points  the  exact  solution  u^x  is 
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FIGURE  24.  The  vertex-edge  neighborhood  ^£5  ^(Ai,  Ai2)* 


generally  not  known  explicitly,  i.e.,  neither  the  exact  eigenpairs  nor  the  exact  EFIFs,  VFIFs  are 
known,  therefore  numerical  approximations  must  be  found. 


3.2.2  Determination  of  the  Eigenpairs 


First,  the  modified  Steklov  method  for  determining  the  eigenpairs  in  the  neighborhood  of  the 
singular  point  is  described.  Herein  the  formulation  for  the  elastic  problem  is  provided  for  the  case 
of  the  edge  singularity  described  in  Section  3.2.1.  The  formulation  for  the  heat-transfer  problem  is 
similar  but  much  less  complex,  therefore  not  presented  herein.  Consider  the  domain  shown  in 
Figure  25,  which  represents  a  cross-section  of  an  edge  singularity,  and  where  r,  0  are  the  coordi¬ 
nates  of  a  cylindrical  coordinate  system  located  in  the  singular  edge. 

By  formulating  the  weak  form  over  Q*r,  the  singular  edge  is  excluded  from  the  domain  of 
interest  such  that  the  accuracy  of  the  finite  element  solution  does  not  deteriorate  in  its  vicinity.  On 
the  boundaries  Fj  and  r2  traction-free  or  zero  displacement  boundary  conditions  are  assumed: 


T’  =  0  or  w  =  0  on  r^*,  1 ,  2,  3 

In  Q*r,  Ux,  Uy  and  may  be  represented  as  follows: 


“.c 

/i(e) 

a 

u  =  * 

Uy 

.  =  r  * 

/2(e)  ■ 

.  “z  . 

.  m) . 

(41) 


(42) 
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Material  3 


Material  2 


Material  1 


Singular  Edge 

FIGURE  25.  Cross  section  of  an  edge  singularity  for  the  modified  Steklov  formulation. 


Using  Eq.  (42)  on  r3: 


and  a  similar  condition  on  E/. 


=  (a/R)u 
dr 


Multiplying  the  equilibrium  equation  by 


V  =  {V^Vy,  V^] 


’^e  X  integrating 


using  Green’s  theorem,  and  following  the  steps  presented  in  Ref.  [1],  the  modified  Steklov  weak 
form  is: 


Seek  oceC,  0  ^  u  g  H  x  H 

fi(M,v)-(iV«(M,v)  +  N^.(M,v))  =  a(M«(M,v)  +  M^.(M,v)),  V  ve  x//'(Q*)  (44) 

where  C  is  the  complex  plane  and 


z 

MR(M,h  =  jj[t\Q][E][R]hr  =  R  dQdz, 


N/u,h  =  JJ[r[G][£]([D‘®>]«)],  =  fi  dQdz, 
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and  [D],  [Q]  and  [i?]  are  given  as  follows: 


[D]  = 


4-0  0 

^  d 

0 

-l 

dx 

(-sinO)^ 

0 

0  4-  0 
ay 

0 

0 

0  0  4- 

0 

0 

0 

dz 

d 

.  a 

0 

i-i-  0 

R 

(-sinO)^ 

dy  dx 

0  Ai. 

0 

0 

dz  dy 

A  0  A 

[dz  dxj 

0 

0 

(-sinO)^^ 

(48) 


[Q]  = 


COS0  0  0  sin0  0  1 

0  sin0  0  COS0  1  0 

0  0  1  0  sin  0  COS0 


and  [E]  is  the  material  stiffness  matrix. 


[R]  = 


COS0  0  0 

0  sin  0  0 

0  0  0 
sin  0  cos  0  0 

0  0  sin0 

0  0  COS0 


(49) 


Remark  1.  The  domain  does  not  include  a  singular  edge,  hence  no  special  refinement  of 
the  finite  element  mesh  is  required. 


Remark  2.  The  formulation  of  the  weak  form  was  not  based  on  the  assumption  that  the  mate¬ 
rial  is  isotropic,  and  in  fact  can  be  applied  to  multi-material  anisotropic  interface. 

The  domain  is  divided  into  finite  elements  through  a  meshing  process.  The  polynomial 
basis  and  trial  functions,  } ,  are  defined  on  a  standard  element  in  the  T|,  ^  space  such  that  -1< 
^  <1,  -1  <  T)  <1,  -1<  ^  <1.  The  entries  of  the  unconstrained  stiffness  matrix  corresponding 
to  B{u,v)  are  given  by  (see  Ref.  [1]): 


z  n; 

For  simplicity,  on  Tj  and  r2  traction-free  boundary  conditions  are  assumed.  Consider  first 
n(u,v)  =  N^(u,v)  +  N^.(u,v) .  The  entries  of  the  matrix  [Nj^]  corresponding  to  the  bilinear  form  Nj^  are 
computed  using  Gauss  quadrature: 
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ST  6 

X  (51) 

i  =  1  /  =  1  u  =  1 

where  W,,  are  the  weights  and  T)  are  the  abscissas  of  the  Gauss  quadrature  points,  respec¬ 
tively,  and  Pit  and  are  matrices  given  explicitly  in  Ref.  [1]. 

The  entries  of  [M]^]  are: 

CO-0  ^  ^  ® 

^  I X  ^  Pil{^s->'^l)ElkPkj(^s’'^t)  (52) 

i  =  1 »  =  1  l,k=\ 

Expressions  similar  to  Eq.  (51)  and  Eq.  (52)  exist  for  the  matrices  [Nj^*]  and 

Denoting  the  set  of  all  coefficients  by  {M;ot)>  set  of  coefficients  associated  with  r3  and 

r4  by  {m^},  the  following  eigenproblem  is  obtained: 


([^:]  -  [Nj,]  -  [N^.]){u,J  =  a([M^]  +  [M^.]){«r}  =  a[M]{u^}  (53) 

The  vector  which  represents  the  total  number  of  nodal  values  in  can  be  divided  into  two 
vectors  such  that  one  contains  the  coefficients  {m^},  the  other  contains  the  remaining  coefficients: 
{utot)^  =  {{^r}^^  By  eliminating  {m,„},  the  reduced  eigenproblem  is  obtained: 

[Ks]{u„}  =  a[M]{u^}  (54) 

Solution  of  the  eigenproblem  given  by  Eq.  (54)  yields  approximations  for  eigenpairs  with  high 
accuracy,  efficiency  and  robustness. 

3.2.3  Extraction  of  the  GSIFs. 

Once  the  eigenpairs  have  been  computed,  they  are  used  for  extracting  the  GSIFs  from  the 
finite  element  solution.  The  procedure  is  as  follows:  First  the  elastostatic  problem  is  solved  over 
the  entire  domain  Q.  by  means  of  the  finite  element  method  based  on  the  displacement  formula¬ 
tion,  thus  obtaining  ufe.  Second,  a  small  sub-domain  around  the  singular  edge  is  considered. 
Define  5/j  as  the  set  of  interior  points  of  a  circle  of  radius  R,  centered  on  the  point  R  Qr  is  defined 
by  Q.n  Rj^x  1^,  where  is  a  segment  along  the  edge:  /^  =  {z  I  zj  <  z  <  Z2)  and  Fr  is  the  circular 
part  of  its  boundary.  See  Figure  26. 

The  complementary  variational  principle  over  Qr  can  be  stated  as: 

Seek  <3q  e  such  that 
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FIGURE  26.  Typical  cross-section  of  an  edge  singularity. 


B^(ao,a,)  =  F^(oi)  V  Oi  €  (55) 

Ec(0.r)  being  the  statically  admissible  space  (see  a  detailed  definition  in  Ref.  [2]),  and  and 
are  given  by: 


S,(ao,cfi)  =  III  OolF]  ^OidQ.dz  (56) 

11  u^[Q]aidsdz  (57) 

where  [£]  is  the  material  matrix,  is  that  part  of  the  boundary  where  the  displacement  vector 

u  is  prescribed,  and  [Q]  is  given  in  Eq.  (49). 

For  the  complementary  weak  form  the  trial  and  test  spaces  are  chosen  to  be  linear  combina¬ 
tions  of  the  eigenstresses,  which  are  computed  from  the  eigenpairs,  using  the  stress-strain  rela¬ 
tionship  and  Hooke’s  law.  In  the  three-dimensional  case  it  has  to  be  recognized  that  the  stress 
intensity  factor  is  a  function  of  z.  Therefore  the  stress  intensity  is  discretized  using  polynomial 
basis  functions  which  are  energy  orthogonal  on  the  interval 
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The  ij-th  term  of  the  compliance  matrix  which  corresponds  to  the  bilinear  form  in  Eq.  (56)  is 
given  by: 


=  11  J  r  Yj  (DuiE^iiDu)‘idrdedz 

,0  0,  l,k  =  l 


(58) 


being  the  real  part  of  the  eigenpair,  and  (Du)’  is  a  vector  corresponding  to  the  j-th  eigenfunc¬ 
tion  and  is  given  in  Ref.  [2]. 

The  eigenstress  tensor,  having  been  derived  from  the  eigenpairs,  automatically  satisfies  the 
boundary  conditions  on  all  boundaries  except  Fr,  so  that  the  linear  form  in  Eq.  (57)  degenerates 
into  an  integral  over  the  cylindrical  boundary  FrX/,  alone. 


Replacing  the  vector  u  in  Eq.  (57)  with  the  approximated  finite  element  solution  on  r3  x  /^, 
upE,  the  7-th  term  of  the  load  vector  corresponding  to  the  linear  form  Eq.  (57)  becomes: 


iFc)j  =  R 


dOdz 

(r  =  R) 


(59) 


where  [Q]  is  given  in  Eq.  (49). 

Solving  Eq.  (55),  one  obtains  an  approximation  for  the  coefficients  of  the  asymptotic  expan¬ 
sion,  the  GSIFs.  Numerical  tests  demonstrated  that  the  rate  of  convergence  of  the  GSIFs  is  as  fast 
as  the  convergence  of  the  strain  energy,  therefore  the  method  is  “superconvergent”.  Analogous 
procedures  will  be  needed  for  the  computation  of  eigenpairs  and  stress  intensity  factors  at  vertex 
and  vertex-edge  singularities. 

3.2.4  Thermal  Loading:  Computation  of  theTSIFs. 

In  [3]  a  mathematical  analysis  on  the  influence  of  the  temperature  field  on  the  stress  intensity 
factors  is  presented.  The  main  results  show  that  elastic  bodies  containing  singular  points  subjected 
to  steady-state  heat  distribution  experience  stress  intensification  determined  by  thermal  stress 
intensity  factors  (TSIFs).  These  TSIFs  are  the  first  two  coefficients  of  the  asymptotic  expansion  in 
Eq.  (1).  The  principle  of  complementary  energy,  together  with  the  modified  Steklov  method  and 
the  p-version  of  the  finite  element  method  can  be  utilized,  as  shown  in  Section  2.3,  for  obtaining 
the  TSIFs  in  a  limit  process  as  the  integration  radius  around  the  singular  point  approaches  zero. 
Although  very  good  results  can  be  obtained,  it  is  necessary  to  refine  the  finite  element  mesh  in  the 
vicinity  of  the  singular  point.  This  is  because  the  present  method  relies  on  the  displacement  field 
close  to  it.  Importantly,  the  proposed  method  is  applicable  not  only  to  singularities  associated  with 
crack  tips,  but  also  to  multi-material  interfaces  and  non-homogeneous  materials.  It  has  been 
shown  also  that  for  weak  stress  singularities  the  method  has  to  be  modified  to  include  terms  asso¬ 
ciated  with  the  displacement  field  due  to  the  singular  temperature  field  around  the  singular  point. 
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Once  these  terms  are  included,  the  integration  radius  can  be  enlarged  to  obtain  good  approxima¬ 
tion  for  the  TSIFs  without  the  need  for  mesh  refinement.  The  implications  of  the  analysis  on  the 
computation  of  the  TSIFs  are  the  following: 

(a)  The  most  singular  term  in  the  stress  tensor  in  the  neighborhood  of  a  singular  point  is  inde¬ 
pendent  of  the  asymptotic  expansion  of  temperature. 

(b)  The  TSIFs  are  independent  of  the  thermal  field  in  the  vicinity  of  the  singular  point.  Never¬ 
theless,  these  TSIFs  depend  on  the  far  thermal  loading,  which  is  the  solution  of  the  linear 
stationary  heat  transfer  problem  affected  by  the  singular  point. 

(c)  The  TSIFs  may  be  extracted  using  the  modified  Steklov  method  and  the  complementary 
energy  method  in  a  limit  process  as  the  integration  radius  approaches  zero,  without 
addressing  the  thermal  distribution  in  the  neighborhood  of  the  singular  point. 

3.3  Schedule 

It  is  estimated  that  the  activities  leading  to  the  completion  of  the  three-dimensional  implemen¬ 
tation  will  take  two  years.  A  Phase  II  STTR  proposal  was  submitted  to  the  Air  Force  Office  of 
Scientific  Research  on  September  12, 1996  with  a  detailed  performance  schedule. 

4  Summary 


All  of  the  objectives  set  for  the  Phase  I  project  have  been  achieved.  Capabilities  for  the  evalu¬ 
ation  of  the  mechanical  strength  of  electronic  components  subjected  to  thermal  and  mechanical 
loading  was  developed  for  two-dimensional  applications.  The  implementation,  based  on  recent 
technological  advances,  makes  it  possible  to  determine  the  natural  straining  modes  and  their 
intensities  at  singular  points  associated  with  multi-material  interfaces. 

The  computational  techniques  were  implemented  in  a  two-dimensional  setting  as  a  module 
within  the  existing  software  infrastructure  of  Stress  Check,  which  was  developed  by  the  small 
business  concern  (ESRD)  over  the  past  six  years.  Stress  Check  is  based  on  the  p-  and  hp-versions 
of  the  finite  element  method  and  has  several  innovative  capabilities  not  available  in  other  finite 
element  programs.  The  new  capabilities  are  available  for  professional  use  through  Stress  Check. 
The  user-information  is  available  in  Chapter  12  of  the  user’s  manual  of  Stress  Check  [17]. 

The  Phase  I  project  also  investigated  the  feasibility  of  the  implementation  of  the  computation 
of  the  eigenpairs  that  characterize  the  temperature  and  displacement  fields  in  the  vicinity  of  singu¬ 
larities  caused  by  multi-material  interfaces  in  three-dimensions.  Successful  completion  of  the 
activities  to  extend  the  current  implementation  into  three-dimensions,  will  lay  the  groundwork  for 
quantitative  evaluation  of  the  conditions  that  cause  mechanical  failure  in  electronic  devices  and 
composites.  This  technological  development  is  an  essential  prerequisite  to  proper  interpretation  of 
experimental  data  in  much  the  same  way  as  the  ability  to  compute  stress  intensity  factors  in  linear 
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elastic  fracture  mechanics  is  an  essential  prerequisite  to  proper  interpretation  of  experimentally 
obtained  crack  propagation  data.  This  analytical  capability,  coupled  with  reference  data  obtained 
from  simple  experiments,  will  make  it  possible  to  evaluate  design  alternatives  in  the  fields  of  elec¬ 
tronic  component  design  and  composite  materials  technology.  The  capability  will  be  made  avail¬ 
able  for  professional  use  through  the  finite  element  analysis  software  Stress  Check.  The  new 
capability  is  expected  to  be  of  substantial  interest  to  manufacturers  of  electronic  components, 
aerospace  companies,  and  their  suppliers. 

The  project  addressed  the  thermo-mechanical  aspects  of  failure  initiation  in  electronic  compo¬ 
nents  under  the  assumption  that  the  natural  straining  modes  of  the  linear  thermoelastic  problem 
characterize  failure.  It  is  expected  that  for  many  materials  in  the  normal  operating  temperature 
range  it  will  be  possible  to  verify  this  assumption  experimentally. 

5  References 


[1]  Yosibash,  Z.  and  Szabo,  B.  A.  Numerical  Analysis  of  Singularities  in  two-dimensions.  Part  1: 
Computation  of  Eigenpairs.  Int.  Jour.  Numer.  Meth.  Engrg.,  Vol.  38,  pp.  2055-2082  (1995). 

[2]  Yosibash,  Z.  and  Szabo,  B.  A.  Numerical  Analysis  of  Singularities  in  two-dimensions.  Part  2: 
Computation  of  Generalized  Flux/Stress  Intensity  Factors.  Int.  Jour.  Numer.  Meth.  Engrg., 
Vol.  39  (3),  pp.  409-434  (1996). 

[3]  Yosibash,  Z.  Numerical  Thermo-Elastic  Analysis  of  Singularities  in  Two-Dimensions.  Inter¬ 
national  Journal  of  Fracture,  Vol.  74,  pp.  341-361  (1996). 

[4]  Babuska,  I.  and  Miller,  A.  The  Post-processing  Approach  in  Finite  Element  Method  -  Part  2: 
The  Calculation  of  Stress  Intensity  Factors.  Int  Jour.  Num.  Meth.  Engrg.,  Vol.  20,  pp.  1111- 
1129,  (1984). 

[5]  Oh,  H.  S.  and  Babuska,  I.  P- Version  of  the  Finite  Element  Method  for  the  Elliptic  Boundary 
Value  Problem  with  Interfaces.  Computer  Meth.  Appl.  Mech.  Engrg.,  Vol  97  (2),  pp.  211-231 
(1992). 

[6]  Szabo,  B.  A.  and  Babuska,  I.  Finite  Element  Analysis.  John  Wiley  and  Sons,  Inc.,  1991. 

[7]  Suo,  Z.  Mechanics  of  Interface  Fracture.  PhD  Thesis,  Harvard  University,  Cambridge,  Massa¬ 
chusetts,  1989. 

[8]  Chen,  D.  H.  Analysis  of  Singular  Stress  Field  Around  the  Inclusion  Comer  Tip.  Engineering 
Fracture  Mechanics,  Vol.  49  (4),  pp.  553-546  (1994). 


40 


References 


[9]  Yosibash,  Z.  Numerical  Analysis  of  Singularities  and  First  Derivatives  for  Elliptic  Boundary 
Value  Problems  in  Two-Dimensions.  D.  Sc.  dissertation,  Sever  Institute  of  Technology,  Wash¬ 
ington  University,  St.  Louis,  Missouri,  1994. 

[10]  Lee,  K.  and  Cho,  Y.  H.  Boundary  Element  Analysis  of  Thermal  Stress  Intensity  Factors  for 
Cusp  Cracks.  Engineering  Fracture  Mechanics,  Vol.  37  (4),  pp.  787-798  (1990). 

[11]  Erased,  N.  N.  V.,  Aliabadi,  M.  H.,  and  Rooke,  D.  R  The  Dual  Boundary  Element  Method  for 
Thermoelastic  Crack  Problems.  International  Journal  of  Fracture,  Vol.  66,  pp.  255-272 
(1994). 

[12]  Sumi,  N.  and  Katayama,  T.  Thermal  Stress  Singularities  at  tips  of  Griffith  Crack  in  a  Finite 
Rectangular  Plate.  Nuclear  Eng.  Desgn.,  Vol.  60,  pp.  389-394  (1980). 

[13]  Anderson,  B.,  Falk,  U.,  Babuska,  I.  and  Von-Petersdorff,  T.  Reliable  Stress  and  Fracture 
Mechanics  Analysis  of  Complex  Components  Using  a  h-p  Version  of  FEM.  Int.  Jour.  Numer. 
Meth.  Engrg.,  Vol.  38,  pp.  2135-2163  (1995). 

[14]  Babuska,  L,  Von-Petersdorff,  T.  and  Anderson,  B.  Numerical  Treatment  of  Vertex  Singulari¬ 
ties  and  Intensity  Factors  for  Mixed  Boundary  Value  Problems  for  the  Laplace  Equation  in 
R^.  SIAM.  Jour.  Numer.  Analysis,  Vol.  31(5),  pp.  1265-1288  (1994). 

[15]  Dauge,  M.  Elliptic  Boundary  Value  Problems  in  Corner  Domains  -  Smoothness  and  Asymp¬ 
totics  of  Solutions.  Lecture  Notes  in  Mathematics  1341,  Springer- Verlag,  Heidelberg,  1988. 

[16]  Guo,  B.and  Oh,  H.S.  The  Method  of  Auxiliary  Mapping  for  the  Finite  Element  Solutions  of 
Elliptic  Partial  Differential  Equations  on  Nonsmooth  Domains  in  R^.  Preprint  submitted  for 
pubhcation  to  Mathematics  of  Computations,  1996. 

[17]  Stress  Check  User’s  Manual,  Release  2.2,  October  1996.  Engineering  Software  Research  and 
Development,  St.  Louis,  Missouri. 


